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ABSTRACl 

Over  the  last  five  years.  so\craI  new  fast  algonthiTr;  lin\c  been  tle\cl(iped  fcir  Ihc  sdlution 
of  elliptic  finite  elemctil  prcihiems  with  many  dctircc^  ef  freedom  Among  them  are 
Yserentant's  hicrarehica!  basis  muitigrid  method  aiui  n  number  ol  domain  ileeomposition 
algorithms  for  problems  ilnided  into  many  subprobleins  I  he  eondition  number  of  the 
relevant  operators,  for  the  best  of  these  prcconditionetl  conjugate  grailient  methods. 
gi^ows  only  slowly  with  the  si/e  of  the  problems,  the  growth  typically  is  quadratic  in  the 
logarithm  of  the  number  of  tlcgrees  of  freedom  associated  with  n  subregion  or  an  clement 
of  the  coarsest  mesh. 

Important  components  of  many  domain  decomposiii(in  precondilioners  are  subproblems 
related  to  the  si-ts  of  s.inables  on  the  interfaces  b^•l^uTn  neigliborine  subreizions  In  this 
paper,  we  deinonstrali'  ih.il  for  problems  in  tuo  iluncn'-ions.  a  very  simple  change  of 
basis  for  these  one-ditncnsional  problems  leads  to  .1  prccondilioiur  uhich  is  ,is  elfcclive 
as  those  prc\'iously  consKkTcd  In  our  proof  we  use  only  tools  of  linear  ilyehra  and  a 
result  of  Ysercntanl's  Our  numerical  expenmenis  ((infirm  lli.at  ihc  new  .ilgonthm  is 
better  conditioned  than  ihc  original  hierarchical  basis  tnulligrul  method 
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1.  Introduction.  We  consider  second  order,  self  adjoint,  uniformly  elliptic  dif- 
ferenticJ  equations  on  a  two  dimensional  polygonal  domciin  Q.  The  problems  are  solved 
numerically  by  using  continuous,  piecewise  linear  finite  elements.  The  domain  is  first 
subdivided  into  nonoverlapping,  triangular  subregions  Q,,  also  called  substructures. 
and  these  are  further  triangulated  into  elements.  H  denotes  the  diameter  of  a  typical 
substructure  and  h  the  diameter  of  its  elements. 

We  develop  a  domain  decomposition  algorithm  similar  to  those  considered  by 
Bjorstad  and  Widlund  [3], [4],  Bramble,  Pasciak.  and  Schatz  [5], [6],  Dryja  and  Widlund 
[13]  and  Widlund  [19],  When  using  these  methods  the  variables  interior  to  individual 
substructures  are  first  eliminated.  The  resulting  reduced  system,  the  Schur  comple- 
ment, therefore  only  involves  the  variables  aissociated  with  F,  the  set  of  edges  and 
vertices  of  the  substructures.  This  system  is  then  solved  by  a  preconditioned  conju- 
gate gradient  method,  where  the  preconditioner  is  constructed  from  certain  problems 
associated  with  the  interfaces  r,j  =  <9Q,  fl  dQj  between  the  substructures,  and  a  global 
coaxse  problem  associated  with  the  vertices  of  the  substructures.  We  note  that  it  is 
shown  in  Dryja  and  Widlund  [13],  that  such  a  preconditioner  naturally  can  be  viewed 
in  terms  of  a  splitting;  cf.  Varga  [18].  In  the  splitting,  the  couplings  between  the 
groups  of  variables,  associated  with  individual  edges  of  the  substructures,  are  elimi- 
nated. In  [13],  we  also  explaiin  how  results  and  algorithms  for  the  two  substructure 
case  cam  be  used  to  construct  and  cmalyze  problems  on  many  substructures. 

Various  preconditioners  have  been  proposed  for  the  subproblems  associated  with 
the  edges  r,_,.  For  each  F.^,  this  is  essentiadly  a  two  subregion  problem,  and  we  can 
therefore  take  advantage  of  a  number  of  results  obtained  in  early  work  on  domain 
decomposition  algorithms.  Already  in  1980,  Dryja  [11],  see  also  [12],  introduced  an 
effective  preconditioner  J,  which  is  the  square  root  of  a  discrete,  one  dimensional 
Laplacizin  on  r,j.  The  same  preconditioner  was  also  discussed  in  Bjorstad  and  Wid- 
lund [4]  and  Bramble,  Pascieik  amd  Schatz  [5].  Other  preconditioners  for  these  two 
subregion  subproblems,  such  as  the  Neumann- Dirchlet  algorithm,  were  considered  by 
Bjorstad  and  Widlund  [4],  Bramble,  Pasciak  and  Schatz  [6],  Chan  and  Resasco  [S]. 
Chain  and  Keyes  [7],  Dihn,  Glowinski,  Periaux  [10]  and  Golub  and  Mayers  [14].  .\ 
number  of  the  resulting  adgorithms  for  the  majiy  substructure  case  are  known  to  be  al- 
most optimad  in  the  sense  that  the  condition  number  is  bounded  by  C{\  +\og(H/h  i  iv 
For  a  more  complete  discussion,  see  Bjorstad  and  Widlund  [4]  and  Widlund  [19]  for 
the  two  subregion  amd  mainy  subregion  caises,  respectively. 

An  adternative  almost  optimal  algorithm,  which  uses  a  hierarchical  basis,  has  been 
introduced  by  Yserentant  [21].  His  bound  on  the  condition  number  is  of  the  same  form. 
When  the  standaird  finite  element  nodad  basis  is  replaced  by  a  hierarchical  basis,  the 
transformed  matrix  becomes  much  better  conditioned.  The  preconditioner  is  then 
obtained  by  discarding  the  off  diagonal  blocks  and  by  replacing  all  but  one  of  the 
diagonal  blocks  by  diagonal  matrices. 

In  this  paper,  we  consider  a  hybrid  method  demonstrating  that  a  successful  and 
simple  preconditioner  can  be  obtained  by  changing  the  bases  of  the  spaces  associated 
with  individual  edges  F,j.  Our  proof  uses  only  tools  of  lineair  algebra  and  Yserentant 's 
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Abstract.  Over  the  last  five  years,  several  new  fast  algorithms  have  been  developed  for  the 
solution  of  elliptic  finite  element  problems  with  many  degrees  of  freedom.  Among  them  are  Yserentant  s 
hierarchical  basis  multigrid  method  and  a  number  of  domain  decomposition  algorithms  for  problems 
divided  into  many  subproblems.  The  condition  number  of  the  relevant  operators,  for  the  best  of  these 
preconditioned  conjugate  gradient  methods,  grows  only  slowly  with  the  size  of  the  problems:  the  growth 
typically  is  quadratic  m  the  logarithm  of  the  number  of  degrees  of  freedom  associated  with  a  subregion 
or  an  element  of  the  coarsest  mesh. 

Important  components  of  many  domain  decomposition  preconditioners  are  subproblems  related  to 
the  sets  of  variables  on  the  interfaces  between  neighboring  subregions.  In  this  paper,  we  demonstrate 
that,  for  problems  in  two  dimensions,  a  very  simple  change  of  bsisis  for  these  one-dimensional  problems 
leads  to  a  preconditioner  which  is  as  effective  as  those  previously  considered.  In  our  proof,  we  use  onl> 
tools  of  linear  algebra  and  a  result  of  Yserentant  s  Our  numerical  experiments  confirm  that  the  n^^w 
algorithm  is  better  conditioned  than  the  original  hierarchical  basis  multigrid  method 
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main  result.  We  show  that  the  new  method  hcis  a  smaller  condition  number  than 
Yserentant's  origineJ  method:  thus  it  grows  no  faster  than  C(l  +  \og{  H / h))'^ ,  a  result 
confirmed  in  our  numerical  experiments. 

Our  work  has  been  inspired  by  recent  work  of  Babuska,  Craig,  Mandel  and 
Pitkaranta  [1],  Babuska.  Griebel,  and  Pitkaranta  [2]  and  Mandel  [17],  [15],  [16]  where 
efficient  preconditioners  for  the  p-version  of  the  finite  element  method,  are  devel- 
oped by  using  hierairchical  basis  functions  and  partial  orthogonalization  of  the  basis 
functions.  Similarly,  our  algorithm  for  the  h-version  involves  a  change  of  basis.  It  is 
extremely  easy  to  carry  out,  and  it  results  in  a  much  better  conditioned  linear  system. 

2.   Yserentant's   Algorithm  and  Iterative   Substructuring  Methods.  In 

this  section  we  provide  the  necessary  background  to  make  it  possible  to  define  the 
new  algorithm  and  to  put  it  in  some  perspective. 

We  consider  a  second  order,  self  adjoint,  uniformly  elliptic,  bilinear  form  aQ(u.v) 
on  Q.  and,  for  simplicity,  impose  a  homogeneous  Dirichlet  condition  on  dQ: 

a(u,v)  =  (f,v),      yveH^iQ).      ue/fo(^)- 

For  the  two  levels  of  triangulations  into  substructures  il,  and  elements  introduced 
earlier,  we  assume  shape  regularity  and  that  they  satisfy  the  usued  rules  of  finite 
element  triangulations;  see  e.g.  Ciarlet  [9].  V'^(fi)  and  V^{Q)  are  the  spaces  of 
continuous,  piecewise  linear  functions,  on  the  two  triangulations,  which  vcinish  on  the 
boundary  dQ. 

The  discrete  problem  is  then  of  the  form 

(1)  a{uH,VH)  =  if,VK),    yv.ev'm,     u.ev'iQ). 

2.1.  The  Hierarchical  Basis  Method.  The  hierarchical  ba^is  method  pro- 
vides a  genersd  purpose  preconditioner  for  second  order  elliptic  problems  in  the  plane; 
see  Yserentant  [21], [20].  The  jdgorithm  is  given  in  terms  of  a  set  of  spaces  V^'.i  = 
0,  •  ■  • ,  J,  which  are  successive  refinements  by  a  factor  of  two  of  V^  =  V^ .  V^'  is  the 
set  of  piecewise  linear  finite  element  functions  after  i  levels  of  refinement  from  the 
original  coarse  triajigulation  with  V^>  =  V''.  V^  is  a  direct  sum  of  subspaces 

v^''  =  Vo''  e  V,''  e  •  •  •  e  v;\ 

where  V',''  =  V^'  \  V'*'-' .  In  other  words,  V,''  is  the  set  of  piecewise  linear  functions  in 
V  ' ,  which  axe  zero  at  the  nodal  points  of  the  triangles  of  aJl  coarser  triangulations,  see 
Figure  1  for  a  one-dimensionaJ  case.  For  the  spaces  V^^,  we  choose  a  basis  of  standard 
nodal  functions  of  V'*'  associated  with  the  new  nodes.  The  resulting  basis  for  the 
entire  space  V^  is  much  closer  to  being  orthogonad  in  the  //'  sense,  than  the  standard 
nodaJ  bsisis  functions,  ajid  the  stiffness  matrix  is  therefore  much  better  conditioned. 

Yserentant's  preconditioner  is  a  block  diagonsd  matrix  in  the  new  basis.  The  first 
block  is  defined  by  the  finite  element  model  for  the  subspace  V^  and  the  others  are 
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Fig.  1     Hierarchical  Basis  Functions  m  One  Dimension 

diagonal.  In  matrix  notation  the  resulting  system,  which  is  solved  by  a  conjugate 
gradient  method,  is  of  the  form 

or  equivalently, 

D-'H^KHi  =  g. 

Here  K  is  the  original  stiffness  matrix,  H  represents  the  transformation  from  the 
hierarchicai  to  the  nodal  basis,  and  D  is  block  diagonal  and  obtained  from  H^ K H . 
as  described  above.  In  an  implementation,  there  is  no  need  to  represent  the  stiffness 
matrix  explicitly  in  the  new  basis.  Instead,  we  perform  the  basis  change  on  the  vectors 
as  needed. 

In  [21],  Yserentzint  develops  the  algorithms  needed  for  performing  the  basis  change 
from  hierajchical  basis  to  nodaJ  beisis  amd  back  ajid  demonstrates  that  each  requires 
fewer  than  2n  additions  and  n  divisions  by  2,  where  n  is  the  dimension  of  the  finite 
element  space.  The  following  algorithm  is  valid  for  both  one  and  two  dimensions. 

Algorithm  to  form  x  «—  Hx 

for  A:  =  1  to  number  of  levels 
for  i  on  level  k 

I,  =  I,  +  (x/i,  +  inMI 
next  i 
next  k 

The  integer  arrays  /I  and  12  contain  pointers  to  the  two  neighbors  of  i  which  are  on 
the  next  coaxser  level.  We  can  regard  the  Jilgorithm  as  defining  a  factored  form  of 
the  matrix  H.  The  noded  to  hierzirchicaJ  transformation  is  similar.  We  note  that  the 
coarse  mesh  needs  not  be  uniform;  see  Figure  1.  If  the  refinements  are  not  uniform. 
then  the  weights  in  the  edgorithms  have  to  be  adjusted.  It  is  also  possible  to  continue 
the  refinement  only  in  selected  subregions. 

2.2.  Iterative  Substructuring  Methods.  Iterative  substructuring  algorithms 
use  a  different  splitting  of  the  space  l'**  into  .V  +  1  subspaces; 

v'"  =  v;1_  e  v'o^n,)  e  •  •  •  e  Vo^^s)  • 
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For  each  substructure  Q,,  we  thus  have  a  subspace  Voin,)  =  V^OH^iQ,).  The 
elements  of  VJj'j,^^  are  piecewise,  discrete  harmonic  functions,  i.e.  they  are  orthogonal, 
in  the  sense  of  the  biUnear  form  a(-,  ■),  to  all  the  other  subspaces.  It  is  easy  to  show 
that  an  element  of  V'/^^^  is  uniquely  determined  by  its  values  on  F  =  [JdQ.,. 

In  a  first  step  of  many  substructuring  algorithms,  the  variables  interior  to  the 
Cl,  are  eliminated.  We  partition  the  vector  x  =  (x/,xs)  and  the  stiffness  matrix 
A'  accordingly.  The  system  that  remains  to  be  solved  is,  after  a  block  Gaussian 
elimination  step, 

(2)  Schur{K)xB  =  9- 

Here  Schur{K)  is  a  Schur  complement  defined  by 

Schur(K)  =  K'bb  -  ^Jb^JI  ^IB- 

A  particular  iterative  substructuring  method  is  defined  by  the  choice  of  a  pre- 
conditioner  for  equation  (2).  Finally,  when  accurate  enough  values  on  F  have  been 
computed,  the  values  elsewhere  are  determined  by  solving  .V  separate  Dirichlet  prob- 
lems on  the  individual  substructures.  We  note  that  it  is  not  necesseu-y  to  compute  the 
elements  of  Schur{K)  since,  in  the  conjugate  gradient  iteration,  this  matrix  is  needed 
only  in  terms  of  matrix- vector  products.  Such  a  product  can  be  found  at  the  expense 
of  solving  one  problem  on  each  of  the  substructures. 

A  number  of  preconditioners  can  be  described  as  follows:  We  first  carry  out 
a  partial  chaxige  of  basis,  associating  the  stajidard  basis  functions  of  I'"  with  the 
vertices  of  the  substructures.  In  the  new  basis,  we  represent  the  Schur  complement 
as 

I  See    Sev 

\   Sgy     Svv 

Here  Svv  denotes  the  part  of  the  Schur  complement  associated  with  the  vertices  of 
the  substructures  and  See  the  past  associated  with  the  edges  between  substructures, 
etc. 

The  preconditioner  for  this  system  is  given  by 

See      0 
0      Svv 

where  Svv  is  the  coarse  mesh  finite  element  problem  atnd  See  is  a  block  diagonal 
matrix.  Each  of  its  blocks  is  associated  with  the  variables  of  a  single  edge  F,j.  The 
operator  J,  mentioned  before,  can  be  used  for  this  purpose;  for  other  exeimples  of 
such  algorithms,  see  the  references  given  in  the  third  pairagraph  of  Section  1. 

3.  The  Hybrid  Algorithm.  We  could  combine  the  two  main  ideas  of  Section 
2  as  follows:  We  first  represent  the  stiffness  matrix  in  the  hierarchical  basis  and  then 
eliminate  the  interior  variables  of  aJl  the  substructures.    We  proceed  by  solving  the 
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remaining  Schur  complement  system  approximately  without  further  preconditioning. 
Finally  we  use  the  resulting  values  as  boundary  data  for  the  local  problems  on  the 
individual  substructures. 

In  the  new  algorithm,  we  proceed  differently,  but  as  we  will  see.  we  will  obtain  the 
same  approximate  solution  without  the  considerable  expense  of  converting  the  stiffness 
matrix  into  the  hierarchical  bcisis.  In  our  algorithm  we  work  with  the  standard  nodal 
basis  while  eliminating  the  interior  variables,  only  changing  to  the  hierarchical  basis 
on  r,  the  set  of  interfaces  and  vertices.  The  resulting  linear  system  is  similar  to  that 
of  Section  2.2 

DgBHggSchur(K)HBBiB  =  g 

It  is  important  to  note  that  we  do  not  use  any  further  preconditioning  of  the  variables 
associated  with  the  edges  F,;. 

This  method  offers  severed  possible  advantages  over  standard  hierarchical  basis. 
The  conjugate  gradient  iteration  is  carried  out  over  a  much  smaller  set  of  unknowns 
and  we  will  show  that  the  condition  number  is  smaller.  The  solution  of  the  subprob- 
lems  is  easily  paxallelizable  since  they  are  independent.  The  hierarchical  bases  method 
in  its  original  form  appears  to  offer  less  opportunity  for  this  trivial  type  of  paralleliza- 
tion.  The  chemge  of  basis  required  in  each  iteration  step  now  consists  of  completely 
independent  one  dimensional  problems  instead  of  a  two  dimensional  problem.  The 
basic  observation  is  that  the  values  at  a  node  on  F.^  can  be  computed  using  only  the 
coefficients  for  the  hierarchical  basis  functions  related  to  that  edge. 

We  now  prove  the  sdmost  optimality  of  our  algorithm  using  two  simple  lemmas 
and  Yserentant's  result. 

Lemma  1.  Let  G  represent  a  change  of  basis  which  leaves  the  space  of  varxahlt's 
on  F  invariant.  Then  the  Schur  complement  associated  with  this  set  of  unknowns  is 
independent  of  the  choice  of  bases  for  Vq{Qi,). 

Proof.  Let  i/  be  the  vector  of  unknowns  associated  with  lQ''(nfc),  V^.  and  xg  be 
those  associated  with  F.  The  most  general  basis  transformation  considered  here  is  of 
the  form 

Gil     Gib    ]  (  ^i 
0      Gbb  J  \  is 

In  the  new  basis,  the  stiffness  matrix  is 

j>^[  GJ,        0     \  (  Kn     K,B   \  (  Gu     G,g 


Gib    Ggg  )  \  Kjg    R'bb  J  \     0      Gbb 
A  straightforward  calculation  shows  that  its  Schur  complement  satisfies 

Schur{K)  =  G^gSchur(K)GBB  ■ 

The  following  result  follows  easily  by  a  Rayleigh  quotient  argument. 
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Refinement  levels 

1 

2 

3 

4 

5 

6 

7 

#  of  unknowns  on  F 

3 

7 

15 

31 

63 

127 

255 

New  method 

1,68 

2.66 

3.82 

5.18 

6.75 

8.52 

10.50 

No  preconditionmg        3.05     6. 88     14.20     28.63     57.37     114.79     230.49 

Table  1 
Condition  Number  for  Two  Subdomam  Case 


Lemma  2.  Let  K  he  symmetric,  positive  definite.  Then,  the  condition  num.beT3 
of  K  and  its  Schur  com.iplem.ent  satisfy 

K{Schur{K))  <  k{K). 

Our  main  result  is  given  in 

Theorem  1.  The  condition  num.ber  of  the  hybrid  algorithm,  introduced  m  this 
section,  is  bounded  by  that  of  Yserentant's  m.ethod.  Thus,  it  is  bounded  by  C(l  + 
\og{H/h))\ 

Proof  We  use  Lemma  1  twice  and  the  fact  that  D  is  block  diagonal  to  obtain 

Schur{D-"^H^KHD-"^)    =    DB'^^Schur{H'^KH)DB'J^ 

=    D-s'J'HlsSchur{K)HBBDs'J'  . 

By  using  Lemma  2,  we  obtain 

K{Dg'J'HlBSchur{K)HBBDBB'}    =    KiSchuriD-'I^HT^  KHD-"^)) 

<    KiD-'I'^H'^KHD-'"'), 

which  is  bounded  by  C(l  -I-  log(if//i))^;  see  Yserentant  (Thm.  4.1)  [21]. 


4.  Numerical  Experiments.  In  a  first  set  of  experiments,  we  consider  the  do- 
main Ct  =  QiUf^2  where  Q,\  zind  CI2  Eire  unit  squares  aligned  eilong  an  edge  F  = 
^1  n  f^2-  We  use  the  standard  regular  mesh  ajid  the  usual  five  point  discretization  for 
the  Laplacian.  The  results  axe  listed  in  Table  1. 

Remark.  Our  experiments  show  that  the  condition  number  grows  faster  than  ( 1  + 
log  H/h)  for  the  two  subdomain  case.  We  note  that  for  a  number  of  preconditioners 
the  condition  number  remains  bounded  in  this  case.  This  is  true  for  the  preconditioner 
J  if  we  solve  a  Dirichlet  problem,  cf.  [4],  but  not  for  a  Neumann  problem.  Yet,  our 
method  and  that  based  on  the  J  operator  both  have  condition  numbers  which  grow 
like  (1  +  log  H/hy  in  the  mziny  subdomadn  case. 

In  a  second  set  of  experiments,  we  consider  the  case  of  majiy  substructures.  The 
unit  squaire  fJ  is  subdivided  uniformly  into  4,  16,  64  or  256  squaire  subdomains  and 
the  same  model  problem  is  solved  using  uniform  meshes.  We  compare  our  results 
with  a  set  of  experiments  reported  in  Yserentant  [21].  His  coarse  problem  hcis  only  a 
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Levels 

3 

4 

5 

6 

7 

S 

9 

#  of  unknowns  in  fl 

r 

15^ 

31^ 

63' 

127- 

255- 

511- 

Yserentant's  method     10.59 
New  method  3.66 

No  preconditioning         9.77 


4  Subdomains 

19.53     31. So  47.14 

6.13       9.37  13.41 

21.50     44.97  91.98 


65.38     86.51     110.49 


#  of  unknowns  on  F        13 


29 


61 


125 


New  method 

No  preconditioning 


#  of  unknowns  on  T 


16  Subdomains 
4.85       7.94      11.80     16.44 
35  75         155        316 


81 


i: 


369 


'53 


New  method 

No  preconditioning 


#  of  unknowns  on  T 


#  of  unknowns  on  F 


64  Subdomains 

5.29       8.52 
137        290 


12.54     17.32 
599       1217 


385        833       1729      3521 


256  Subdomains 

New  method 

5.39 

8.70 

12.78 

17.60 

No  preconditioning 

546 

1152 

2372 

4766 

1665      3585      7425      15105 


Table  2 
Condition  Numbers  for  the  \{any  Subdomam  Case 


single  node  and  it  is  therefore  directly  comparable  only  with  our  experiments  on  four 
substructures.  The  results  are  given  in  Table  2. 

The  coarse  problem  and  the  problems  associated  with  the  edges,  which  together 
make  up  the  preconditioner,  are  independent.  We  can  therefore  scale  the  contribution 
of  the  coarse  model  by  a  scalar  factor  a  selecting  the  \-alue  of  the  parameter  for  which 
the  convergence  is  fastest.  In  our  numericai  experiments,  we  have  found  that  for  our 
model  problem  q  =  3.6  is  the  best  for  a  wide  range  of  refinements.  In  the  case  of  4 
subdomains  with  a  coaxse  problem  consisting  of  a  single  point,  we  have  found  that  the 
best  a  varied  considerably  with  the  mesh,  ranging  from  5.4  to  7.S9.  We  note  that  the 
condition  number  grows  quadratically  in  the  logarithmic  factor  for  any  fixed  a  >  0. 
Our  numerical  results  are  reported  for  optimal  choices  of  q. 

In  other  experiments,  which  will  be  reported  elsewhere,  we  have  found  that  there 
is  very  little  difference  in  the  performance  of  the  new  and  several  previously  known 
domain  decomposition  methods.  We  also  note  that  a  version  of  our  algorithm  has 
been  tested  successfully  for  a  membrane  model  of  plane  linear  elasticity  by  Anders 
Hvidsten  and  Petter  Bjorstad  of  the  University  of  Bergen.  Norway. 
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